
function [q,Dq] = Q_cfs(gamma,ES0,ES1,ES2,W,n)

N0 = n(1,1); N1 = n(1,2); N2 = n(1,3);

g = [ES0(1:N0,:);ES1(N0+(1:N1),:);ES2(N0+N1+(1:N2),:)];
g = [mean(g-ES0,1);mean(g-ES1,1);mean(g-ES2,1)] + gamma;

q = g' * (W \ g);

if nargout > 1
    Dq = 2 * (W \ g);
end



